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Motivated by the recent investigations on instabilities caused by Schwoebel barriers during 
growth and their effects on growth or sublimation by step flows, we have investigated, using the 
Stillinger- Weber potential, how this step edge barrier arises for the two high symmetry steps 
on lxl reconstructed Si(lll). Relative to a barrier of 0.97 ± 0.07 eV on the surface, we find 
additional (Schwoebel) barriers of 0.61 ± 0.07 eV and 0.16 ± 0.07 eV for adatom migration over 
the [211] and the [112] steps respectively. The adatom potential energy is found to be strongly 
correlated with that derived from the local geometry of atoms on the adatom-free surface or 
step edges. This correlation preserves a strict correspondence between the barrier determining 
features in the spatial variation of the adatom potential energy and the same derived from 
the local geometry for the Si(lll) surface and the [211] step. It is therefore argued that the 
Schwoebel barrier on the [211] step is robust i.e. a feature that would survive in more satisfactory 
ah initio or tight binding calculations. Using a diffusion equation for the adatom concentration 
the relevance of the barrier to electromigration of steps has been explored. Data from such 
experiments on Si(lll) has been used to place an upper bound on the Schwoebel barrier and a 
lower bound on the electromigration force. 



PACS numbers: 68.35Fx, 68.35Ja, 68.35Bs, 68.35Md 
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I. INTRODUCTION 



The Schwoebel barrier was originally introduced in the context of step motion ffl as the 
additional barrier for adatom diffusion over a step edge from the upper to lower terraces. 
It was argued that such a barrier results in an anisotropy in adatom diffusion into the 
step edge - the diffusion from the lower terrace being greater. This anisotropy was found 
to drive an arbitrary distribution of step spacings towards a uniform distribution during 
growth of a vicinal surface. JTJ] Later, it was pointed out by Villain P| that this growth by 
step flow is stable only on a sufficiently vicinal surface with possible instabilities setting 
in during the growth of a flat (singular) surface. The dynamical morphology of a singular 
surface growing under the influence of a schwoebel barrier is a subject of great current 
activity. While the eventual fate of growth on a flat surface in the presence of such 

extra step edge barriers is still being discussed , it is now well accepted that Schwoebel 
barriers lead to coarsening in the evolving surface morphology under nonequilibrium growth 
conditions, producing mounds, pyramids and facet-like angular structures on the growing 
surface. Recent experimental studies of nonequilibrium growth on the Si(lll) surface have 
produced somewhat contradictory p|-|10|j results, and the specific role of Schwoebel barriers 
for nonequilibrium Si(lll) growth is unclear at the present time. 

The current study is primarily motivated by observations of another instability: As was 
pointed out by Schwoebel and Shipsey Q, an anisotropy favoring diffusion into the step 
edge from the upper terrace (possibly due to larger barriers for diffusion from the lower 
terrace) results in a step pairing instability during the growth of a vicinal surface. In recent 
experiments similar direct-current induced reversible step-bunching instabilities have been 



observed during sublimation of the high temperature lxl phase of Si(lll). |TT],|T2[- Here, 
we report on a calculation of the Schwoebel barrier for the two high symmetry vicinal steps 
on lxl Si(lll) using the empirical Stillinger- Weber potential. |13| The importance of the 
Schwoebel barrier in the context of the electromigration experiments has also been explored 
by modeling the barrier as affecting the boundary conditions to the diffusion equation for 
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adatom concentration used to interpret these experiments 



II. EMPIRICAL POTENTIAL CALCULATIONS 



The use of the Stillinger- Weber potential in this study has been motivated by its suc- 



cessful application in previous studies of bulk and liquid silicon ||I3] , the Si(100) surface and 
steps on this surface ||15|| . Barriers on the Si(100) surface and on single and double height 
steps on this surface have also been calculated using this potential. [[L8]-[2(J Although this 
potential fails to reproduce the correct energetics of the Si(lll) surface configurations with 
adatoms, Jl6| it has been used here since features that follow purely from the changes in 
the adatom coordination number are expected to be robust i.e. these features survive even 
if the details of the empirical potential used change. Such features are expected to survive 
in more satisfactory ab initio or tight binding calculations. The calculation here is followed 
by an attempt to identify such robust features. 

To determine the diffusion barriers the adatom potential energy has been mapped as a 
function of the (x, y) position of the adatom (in the (111) plane) for the Si(lll) surface (Fig. 
2(a)), the [211] step (Fig. 2(d)) and the [TT2] step (Fig. 2(g)). The three-fold and reflection 
symmetry of the Si(lll) as shown in Fig. 1 implies that steps running along directions with 
equal 9 are identical. It has been shown in a previous study |H]] using the Stillinger- Weber 
potential that an alternative configuration of step-edge atoms with some of the upper terrace 
atoms rebonding to atoms in the lower terrace gives a lower step energy for the [112] and 
[101] steps. However, such a configuration has not been considered here since it has also 
been shown that it gives rise to step-step interactions an order of magnitude larger than 
experimental estimates. pl[ Neglecting such rebonded configurations, all the intermediate 



low symmetry steps (0° < 9 < 60°) have been shown to have a higher step energy. |2T| In 



other words, diffusion barriers have been calculated for those steps whose interactions are not 
larger than experimental estimates and whose step energy is a local minimum as a function 
their orientation 9. The adatom potential energy V has been computed as the difference in 
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the minimum potential energy of the system with the adatom at infinity (non-interacting) 
and the same with the interacting adatom. 

Standard molecular dynamics (MD) procedures of integrating Newton's law (with dis- 
sipation to reduce temperature) and the steepest descent equations have been used to de- 
termine the minimum potential energy of the system. [pTJ These routines determined the 



adatom potential energy to an accuracy of 10 -4 eV. The (x, y) coordinates of the adatom 
are fixed during the integration process. The system consisted of a certain number of bi- 
layers of Si(lll) in an MD cell with the surface lattice constants a\ and 02 along its x and 
y directions respectively (see Fig. 1). Three of the bottom bi-layers are fixed at bulk lattice 
coordinates throughout the calculation. In simulations on the Si(lll) surface with six ad- 
ditional (movable) bi-layers it was found that changing the adatom's (x, y) position from a 
deep minimum lead to the shearing of the entire surface in the xy plane towards the adatom. 
As this made the (x, y) position of the adatom relative to the surface ill defined further sim- 
ulations on the surface as well as the step configurations were carried out with atoms at the 
(x, y) boundaries also fixed at positions corresponding to the adatom-free but relaxed con- 
figurations. The regions in the (111) plane explored in all cases were "centrally" located i.e. 
maximally away from the boundaries so that finite size effects are minimized. System size 
dependence in V(x, y) was initially explored for the Si(lll) surface to determine the optimal 
system size. It was found that changing the system size from four lattice constants in the x 
and y directions (including the boundary of fixed atoms) with three movable bi-layers to six 
lattice constants in the x and y directions with six movable bi-layers changed (reduced) the 
adatom potential energy (at the minima, maxima and saddle points) by < 10~ 2 eV. As an 
error bar of ±10~ 2 eV was estimated to be smaller that the accuracy needed for this study 
all the simulations were carried out with the smaller system size - the size in the y direction 
being extended to 4| and 4| lattice constants for the [211] and [TT2] step configurations 
respectively so that periodic boundary conditions could be applied to create vicinal steps. 
Since the bulk terminated Si(lll) surface and [211] step configurations do not relax under 
the Stillinger- Weber potential, with the system sizes chosen here their corresponding surface 



and step energies are reproduced exactly. However, the [112] step configuration is different 
from the bulk terminated structure due to the 2x1 reconstruction at the step edge. With 
the system size chosen here it's step energy is reproduced to within 4xl0 -4 eV/a\. f2"I|| 

The MD procedures began with initial configurations for each (x, y) position of the 
adatom corresponding to the relaxed adatom-free structures. The z coordinate of the adatom 
was chosen to be equal to that obtained in the final configuration during simulation with 
the adatom in the immediate neighborhood of the point {x,y). For the step configurations 
V(x, y) was determined from behind to the front of the step edge. These procedures for 
determining V(x, y) were found necessary for the [211] step configuration since other methods 
such as an arbitrary initial z coordinate for the adatom or determining V(x, y) from the front 
to behind the step edge lead to the adatom relaxing into configurations in which it displaces 
an atom near the step edge and/or moves into the bulk. Symmetries in the (111) plane 
were exploited to reduce the size of the regions in which the adatom potential energy needs 
to be computed. Therefore for the Si(lll) surface only a sixth of the surface unit cell was 
explored. For this case V(x, y) was computed on a triangular grid - the sides of the triangle 
coinciding with the high symmetry directions with the distance between neighboring grid 
points being An interpolation scheme respecting the symmetries on the surface was used 
to determine the features of V(x,y). For the [211] and [112] step configurations reflection 
symmetry about the y axis (_L to the step edge) reduced the width of the region (along the 
step edge) explored to a length of y and a\ respectively. For these configurations V(x, y) was 
determined on a rectangular grid with the distance between neighboring grid points being 
Yjt and |{j along the x and y axis respectively. With the step edge in the middle, the length 
of the region explored was l|ct2 and l|ci2 for the [211] and [112] steps respectively. This was 
to enable a new interpolation scheme (respecting the symmetries of the step configurations) 
applying periodic boundary conditions along the x axis to do the same along the y axis 
of the region explored. On comparing the potential energy of the adatom on the Si(lll) 
surface to that far away from the [211] step edge errors due to the finite grid sizes, the 
interpolation schemes and possibly effects due to being close to the boundary of fixed atoms 
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were recognized. These errors are conservatively estimated to be ±0.05 eV. Since this is 
much larger than the errors due to finite size effects it is assumed to be the error bar in the 
potential energy. Barrier values being differences in these potential energies are therefore 
estimated to have an error bar of ±0.07 eV. 



III. RESULTS AND DISCUSSION 

The results are shown in Fig. 2. The global minimum of the adatom potential energy V 
occurs on the H 3 site (mi) on the Si(lll) surface (Fig. 2(b)) where the adatom potential 
energy is -3.31 eV. This is a significant fraction of the bulk energy per atom: -4.34 eV. The 
relevant saddle point for H3 <-> H3 transitions is si (close to the T4 cite which is a local 
minimum) where the 

potential energy is -2.34 eV. Thus the barrier to diffusion on the surface is 0.97±0.07 
eV. These results are consistent with previous studies (using the same potential) on surface 



energies of configurations with adatoms [I6| as well as a study of diffusion on the Si(lll) 



surface. [22| The barriers to diffusion between the minima along the step edges are slightly 



smaller than the barrier on the surface - apparently inconsistent with previous work |3~4j 



showing that step edge fluctuations are (predominantly) due to attachment/detachment of 
adatoms from the terrace and not due to diffusion of atoms along the step edge. In this 
study the Schwoebel barrier has been defined to be difference between the maximum adatom 
potential energy (along the path on which this is a minimum) as it moves into the step edge 
from the global minimum on the upper terrace far away from the step edge and the same 
for the H 3 «-> H 3 transitions on the Si(lll) surface. In other words it is the difference in the 
adatom potential energies at the barrier determining saddle point on the step configurations 
and the same on the free surface (si). With this definition the Schwoebel barrier cannot be 
negative. From Fig. 2(e) (Fig. 2(h)), for the [211] ([TT2]) step the barrier is determined 
by S3 (S5) where the adatom potential energy is -1.73 eV (-2.18 eV) - implying a Schwoebel 
barrier of 0.61±0.07 eV (0.16±0.07 eV). Growth on Si(lll) is therefore expected to produce 
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moundlike structures with facets consisting (predominantly) of [211] steps. However, the 
experiments of Yang et al (temperature = 275±5° C) instead show facets with [112] steps 
||. This discrepancy may be due to the presence of the 7x7 reconstruction. 

It must be noted that the Stillinger- Weber potential has been tuned only to the properties 
of the bulk diamond and liquid structures of Silicon and not to any surface or step properties. 
As mentioned previously it does not produce the correct energetics for the Si(lll) surface 
configurations with adatoms |16| as compared to ab initio calculations. |23|-p5| Therefore an 



attempt has been made to identify the robust features of this study: Features following from 
changes in the coordination number of the adatom. This idea is supported by the observation 
that the reconstruction energy of the Si(100) surface calculated using the Stillinger- Weber 
potential [^6[ (0.85 eV) agrees with ab initio calculations [p7[] (0.84 eV). The Tersoff and 
Dodson empirical potentials also give this energy to be of the same order of magnitude. 
p6| The Stillinger- Weber also reproduces the correct order of energies of the [211] and 
[TT2] steps (per step edge atom these values are 0.72 eV and 0.62 eV respectively [ pT| ) 
as well as the presence of rebonding at the [112] step edge as compared to tight binding 
calculations (per step edge atom these values are 0.70 eV and 0.38 eV respectively [|T7jj). 
During reconstruction, the coordination number of atoms (on the (100) surface and [TT2] 
step edge) changes from two to three. It is therefore expected that features following from 
such a change in coordination number are not artifacts of the empirical potential used. 
Here, the adatom energy recomputed without additional relaxation of other atoms due to 
the presence of the adatom is assumed to be a good measure of the coordination number. 
Although this measure is very similar to the actual adatom energy V(x,y), it helps to 
identify features that follow from the geometry of atoms (locally around to adatom) on the 
relaxed adatom-free surface or step edges. Features that do not change significantly due 
to additional relaxations in the presence of the adatom are expected to be robust. The 
adatom potential energy Vi g (x,y) has therefore been recomputed with other atoms fixed at 
positions corresponding to the relaxed adatom-free structures. These results are shown in 
Fig. 2(c) for the Si(lll) surface, Fig. 2(f) for the [211] step and in Fig. 2(i) for the [TT2] 



step. Similarities in the contour plots of V and Vi g suggest a strong correlation between 
them. This correlation is evident from the linear relationship between Vi g and V (Fig. 3) 
The lines in Fig. 3 are best fits to Vi g (x,y) vs V(x,y) which were explicitly computed at 
the grid points. There is no repetition of points that are equivalent due to the symmetries 
of the relevant configurations. Further, for the step configurations, the length of the region 
along the y axis corresponding to these plots is |a 2 on either side of the step edge. The 
rough configuration independence of the relationship between Vi g and V suggests that it is 
a characteristic of the [111] surface. A similar behavior may be true for the Si (100) surface 
and single height steps on it. In this case, borrowing the values of V corresponding to the 
features (minima and saddle points) from previous work by Roland and Gilmer |18|-p0|] , and 
computing V\ g for the same features here, it is found that V\ g = (1.3 ± 0.2)V + (1.2 ± 0.5) 
for the Si(100) surface and Vi g = (1.1 ± 0.1)V + (0.6 ± 0.3) for the combined data from the 
three single height steps. 



A. Relevance to electromigration experiments 

Recent experiments on the electromigration of steps on Si(lll) |Tl~l , |l2ll can be reinter- 
preted in the presence of a Schwoebel barrier by modifying the boundary conditions to 



the diffusion equation for adatom concentration used previously by Stoyanov et al |1J| to 
describe step bunching (see appendix A). The particular modification is to suppress the 
strength of the step as a source of adatoms onto the upper terrace relative to the same onto 
the lower terrace by a factor e e where e = |^ with E s being the Schwoebel barrier. With 
this modification the equations for adatom concentration on a particular terrace has four 
length scales - the diffusion length A and the scale introduced by the electromigration force / 
(=-^) both of which are parameters entering the diffusion equation and two other scales to 
and t\ entering the boundary condition to this equation at the upper and lower terrace step 
edges respectively (to,i = J?^ , a ~ 2 is the density of atoms in a terrace, D is the diffusion 
constant, n e is the equilibrium adatom density and (3 is the step kinetic coefficient - here the 
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Schwoebel barrier is modeled as making this parameter assume different values at the upper 
and lower terrace step edges). These are in addition to the scale introduced by the terrace 
width W. Since the modification is further designed to keep the total strength of the step 
as a source of adatoms a constant independent of E s , to + t\ — t - a constant independent of 
e for fixed k^T (variation of t with E s is included in the estimated range of t ) with ^ = e e . 
The (linear) diffusion equation is then solved, with n e appearing in the boundary conditions 
as a "source" , to obtain the spatial variation of the adatom concentration. 

Using the above solution the equations for the velocities of an array of steps with positions 
Xi are developed- the step index i increases in the step down direction which is also the 
positive x axis. The time scale that enters these equations is the lifetime of an adatom r^. 
These equations have the following form (see appendix A): 

X n = —X[g-(e, tX, fX, WnX- 1 ) + g + (e, tX, fX, W^X" 1 )} (1) 

where 9 e = n e a 2 and the terrace width W n = X n+ i — X n . From the above, equations 
for the rate of change of the terrace widths can be obtained. Linear stability analysis of 
these equations around an average terrace width w would predict when the step bunching 
instability occurs. The linearized equations may be written as follows: 

5 n = k_5 n+1 + (k + - k_)5 n + k + 5 n _i (2) 

where 5 n = W n — w and k± = ^\[^M>\w=w- The instability occurs if J28| k + — k_ > or 



equivalently (since k + + A;_ < for w > 0), when the anisotropy ratio p — < 1. This 
analysis is similar to the work of Ghez et al p9 which includes an external flux with the 



electromigration force being absent. 

The recent experiment of Williams et al on reversible step bunching on Si(lll), which 
measures an "effective" anisotropy ratio p e ff (near bunched steps) in the temperature range 
1155-1215° C (this range includes the corrections due to the emissivity of the optical pyrom- 
eter used in the experiment 10), can now be used to to determine an upper bound on E s 



{Eg) and a lower bound on F (Fe) by solving the equation r 1 = k + — pk^ = 0. The solution 



to this equation is obtained by estimating the parameters w, t and A (in the temperature 
range of interest: 1155-1215° C) and then determining / as a function of e. The estimates 
of the upper and lower bounds on t and A are made conservatively - the range is made 
wide enough so as to include the "true value" . Experimental observations of Latyshev et al 



TI[ show that around a temperature of 1200° C the step velocity varies linearly with the 



terrace width upto a width of 2 fim . A lower bound on the diffusion length |32| of 0.5 



\im is obtained by reproducing this result using equation 1. This estimate is robust against 
large changes in the parameters t and /. To make this estimate, e is conservatively chosen 
to be zero. Further, to first order in toA -1 , g+ + <?_ = — wX^ 1 - independent of t, e and /. 
Therefore the magnitude of the slopes of the above curves correspond to the evaporation 
rates r = An estimate of the upper bound on A is obtained by estimating the diffu- 



sion constant g2[ D (using D = b 2 ue~ Ea / kbT , 6=3.84,4, i/=10 13 s~ 1 , E a =0.97 eV (calculated 



here)) and an estimate of the upper bound on rj/ = (The temperature dependence of r 
in the data from Latyshev et al |TT] is consistent with the theory of Burton et al |30| (BCF) 



with an activation energy equal to the cohesive energy of silicon (Ef,=4.34 eV)). The upper 
bound on 8 e is assumed to be 0.167 - the primary motivation for this is that it is measured 
to be « 0.1 around 900° C and [31] its activation energy E n , which is equal to Ej, — E T (E T is 



the activation energy for r^ 1 ) from the BCF theory, is estimated to be negative since total 
energy calculations [3j| predict smaller surface energies for silicon surfaces with adatoms as 



compared to the lxl surface thereby making E T > E^. The specific value of 0.167 is chosen 
as it corresponds to the density in any -^3x^3 configuration of adatoms wherein all the 
floating bonds on the lxl substrate are saturated by bonding with the adatoms - at higher 
densities adatom interactions are expected to be significant. This bound is therefore required 
here as a measure of internal consistency in the analysis since the diffusion equations used 
correspond to free adatoms. The upper bound on r// is calculated using this value of 9 e . 
As the variation in n e has been argued to be small E n « and E T « Et > E a . The 
estimate of the upper bound on A (=(Dr(j) 1 / 2 ) is therefore made at the lowest temperature 
of interest and is equal to ~ 70 \im. The upper bound used in this study is ~7 times this 
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value, making 0.5 pm < A < 0.5xl0 3 pm. 

The estimate of the range of t is obtained from the range of A, the measured r(T 2 ) [DJ 
and 9 e (Ti) [^TJ and the relation between (3 S = Po + Pi and the measured step mobility r(Xi) 
p4| (Ti corresponds to 900° C and T 2 corresponds to temperatures of interest). At the same 
temperature, (3 S = ^ (see appendix B). From a model of attachment-detachment at the 
step edges (see appendix B) the temperature and 9 e dependences of f3 s is determined: (3 S oc 

-2E a -2{Eg + E a ) 

^—^ — + \_ e b — • The upper bound on f3 s is derived using the above dependence from the 

lower bound = r 3 (i? a =0.97eV, E s is set equal to a value that gives the largest increase). 
The parameter t is then evaluated using t = ■ s ^-. This gives upper and lower bounds 

r(T+)(l-e e (T!) 



on t- t , - F(Tl) 
on 1. i U)i - a 2 A 2 r(T +) 



V 



r(T-)(l-9 e (T+)) 



/ r l„ i -\ 

1, 



e k b T 2 T +Q e (r+) 



1 - where T 2 (T 2 ) corresponds to 

the highest (lowest) temperature in the range of interest. The values of r(T 2 + ) and r(T 2 ~) 



are obtained from the data of Latyshev et al [11] using £'{,=4.34 eV, T(Ti) from a previous 
measurement of Bartelt et al ||34|| , i? a = 0.97 eV, ^ e (T!)=0.1 (measured by Yang e£ a/ pT[ ) 
and 0(T 2 + )=O.167 (as assumed previously). This results in 77x ff^ m < t < 7 - 7x ff^ m x 10 3 - 1 . 
The range of t therefore depends on the value of A. 

The value of the parameter w used here is equal to the average terrace width in the 
experiment of Williams et al [0] (=0.15 [im) who determine p e ff to be 0.20±0.03. It must 
be noted that the extraction of p from the experiment had used a theory |35] which makes 
p e ff = plwX^ 1 = 0) (w was much smaller than the average terrace width since in the theory 
used it corresponded to the distance between bunched steps). In this study it is assumed 
that p e ff = p(w\^ 1 7^ 0) - this may be a reasonable approximation since here X min —0.5pm 
and the W max in the experimental data used to determine p e ff was as large as 0.6 pm. This 
approximation is motivated by the need to study the solution of r _1 = including the full 
non-linearity in wX^ 1 . 

The procedure for determining E% and Fi assumes that an electromigration force F 
causes the step bunching instability (p=0.2) with — F restoring the stability (p > 1). This 
assumption is motivated by the observation |j| that a current in the step up direction 
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causes the step bunching whereas an equal magnitude in the step down direction results 
in uniformly spaced steps. & Fg are determined by studying the solutions to r _1 = 
(fixing t, A and w) for p=0.2 and p=1.0 . The variation of r~ x with F (with typical values 
of t, A, w, E s & p) is shown schematically in inset (a) of Fig. 4. Of interest is the solution 
that exists even in the limit wX^ 1 = (the other two solutions do not exist in this limit). 
The variation of this solution with E s (fixed t, A and w) for p=0.2 and p=1.0 is shown 
schematically in inset (b) of Fig. 4. The two curves show that beyond a certain value of E s 
a force — F (with F resulting in p=0.2) cannot restore stability. This value is therefore the 
upper bound (E™) on E s . Below E™ the magnitude of F needed for p=0.2 is larger. The 
value of F at E% is therefore the lower bound (Fi) on F. F = F(E S = 0,p = 0.2). 

Fig. 4 shows the variation of Fg, E™ and F with the parameters t and X (w = 0.15pm). 
As noted previously the range of the t depends on the value of A, with additional restrictions 
on the upper bound of t. These restrictions are due to the absence of the solution (of 
interest) to r _1 = for small A and large t. This seems to correspond the behavior observed 
at higher temperatures [15] wherein a current in the step down direction causes the step 



bunching instability. For tw < 1, the curves in Fig. 4 correspond to the limit wX 1 — ► 
and can be obtained by solving r _1 = to zeroth order in wA™ 1 . This gives 

r k b T ( (1 + e- 2E °l k » T ) - p(l + e 2E ^ T ) \ 



XH \ p-1 
Using the above expression, E%, Fg and F are found to be 



(3) 



^ = ^n(A±^ (4) 



,3po + l, 
k b T ( 4(1 - pi) 



XH V(3 + po)(3p + l) / 



(5) 



with po = 0.2 (here). Equation 4 gives .Rj=0.05 eV for the experiment of Williams et al 



1~2]| . This is much smaller than the barrier calculated here for the [211] step (0.61±0.07 eV) 
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suggesting that [112] steps were seen in the experiment. From equn. 6 and the electric field 
in the experiment |I2] (700 V/m) the (maximum) charge needed on the adatoms (in units of 
the (negative) electronic charge) is estimated to be reasonably small- in the range 3 x 10~ 2 
to 3 x 10~ 51 . From Fig. 4 it can be seen that with a Schwoebel barrier the charge needed 
is always smaller than with a zero barrier. For tw > 1, the deviation of the curves from 
the above expressions is significant - this is due to the dependence of r _1 on wX^ 1 which 
require that terms of order u>A -1 and higher be considered. In this range, the minimum 
charge needed on the adatoms (using Ft) continues to have a reasonable lower limit as it is 
in the range 2 x 10 1 to 7 x 10~ 5 . 



IV. CONCLUSIONS 

In summary the adatom potential energy contours have been calculated for the lxl 
reconstructed Si(lll) surface and the two high symmetry single height steps on it using 
the empirical Stillinger- Weber potential. From these plots diffusion barriers along the path 
of minimum barrier height can be calculated for transitions between any two minima. The 
results show that barriers for diffusion in the trench along the step edge are smaller than that 
between that global minima on the free surface. There is also a strong correlation between 
that adatom potential energy and the potential energy derived from the local geometry 
of atoms on the adatom free surface or surface or Si(lll) surface (0.97±0.07 eV) and the 
Schwoebel barrier on the [211] step (0.61±0.07 eV) are robust features due to changes 
in adatom coordination number. Interpreting recent electromigration experiments in the 
presence of a Schwoebel barrier shows that smaller values of electronic charge on the adatoms 
can account for the observations. It also suggests, in the limit where the diffusion length (A) 
is much larger than the terrace width (w) , that [112] steps were observed in the experiment 



12|j . However, this conclusion relies on assumptions regarding the nature of the atomic 
processes at the step edge (that these processes may involve two atoms (at a kink site) 
which see the barriers as computed here) and that wX^ 1 — > 0. Future studies addressing 
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these assumptions may interpret electromigration experiments differently. 
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APPENDIX A: ANALYSIS OF THE DIFFUSION EQUATION FOR STEP 

MOTION 



The diffusion equation for non-interacting adatoms subliming in the presence of an elec- 
tromigration force F perpendicular to the step edge and the absence of an external flux 
may be written in the "adiabatic approximation" (which neglects the time derivative of the 



density) as follows: fl~4],|29 



d 2 n ^dn n 
dx 2 dx A 2 



(7) 



where n is the density of adatoms, / = F/k^T, A is the diffusion length and the positive x axis 
coincides with the step down direction. This equation can be solved using the boundary 
conditions on the current j (=D (— ^ + fnj, where D is the diffusion constant). These 
conditions determine the step velocity V ( =a 2 (—ji + j u ), where a 2 is the inverse terrace 
density, j u is the current towards the step edge from the upper terrace and ji is the current 
away from the step edge on the lower terrace) and define the step kinetic coefficient (3. In 
the presence of a Schwoebel barrier (E s ), this parameter is modeled as assuming the values 
/3q and Pi at the upper {x = 0) and lower (x = W) terrace step edges respectively, with 



1^ = e e , where e = The boundary conditions can now be written as follows: 



dn 

~T + f n 

dx 



te e 



x=0 



l + e e 



(n(0) - n e ), 



dn 

~r + f n 

dx 



J x=W 



l + e e 



(n(W) - n e 



(8) 
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where n e is the equilibrium density and t = ^tP 1 . The solution to the diffusion equation 



under the above boundary conditions is shown below: 



n(x 



) = n e (Ae wx ~ l71+ +Be wx ~ lr '-) 



(9) 



with 

V± = fX/2 ±7(1 + (/A/2) 2 ), A 
where 



#e w - 1 '- W ffA " 1 "+' 



5 



1 _ ( 2± -/A)(l+e^ 



(t?— /A)(l+e') 



(>? + -/A)(l+ e g ) 



and ■?/> 



fa— /A)(l+e') 



iAe e ' V ■*■ iAe f ' ^ ^ tX ""^ ^ " 1 tX 

The velocity (V n ) of a particular step edge (X n ) can now be calculated: 

V n = a 2 (-j { (= i(x = 0, W = X n+1 - X n )) + Ju (= j( x = W,W = X n - X n _ x )) 
From the functional form of the density and using D = — {jif is the lifetime of an adatom) 

If 

it may be seen that 



V, 



n P a 



Kg- + 9a 



(10) 



where g± are dimensionless functions of dimensionless arguments given by: 



g4e,tXJX,WX- 1 )=t\ 



e e {A+B-l) 
l+e e 



g+(e,t\,f\,W\- 1 )=t\ 



' 1 "++Be w/A ~ 1 "--1 



l+e e 



with the symmetry g + (—a, b, — c, d) = p_(a, 6, c, rf). 



APPENDIX B: AN ATOMIC MODEL FOR THE STEP KINETIC COEFFICIENT 



The step kinetic coefficient (3 was introduced by Chernov |36] through the supposition 
that the currents (ji from the lower terrace and j u from the upper terrace) towards a step edge 
(at x = X) are proportional to the deviation of the adatom density n in the neighborhood 
of the step from the equilibrium density n e . From the model adopted here it will be shown 
that in the presence of a Schwoebel barrier this parameter assumes different values i.e. /3q 
and (3i determining ji and j u respectively: 



A) 



0(x + ) - e e 



Ju 



0i 



e(x-) - e e 



(ii) 
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where a 2 = l/(terrace density), # = na 2 and 9 e = n e a 2 with the positive x axis in the step 
down direction. These equations serve as the boundary conditions used in this work (eqn 
8). 

The temperature dependence of Pq (Pi) can be derived assuming an activated model for 
adatom attachment-detachment at the step edge from the lower (upper) terrace with the 
single atom activation barrier being E a (E a + E s ) where E a is the barrier on the free surface 



and E s is the Schwoebel barrier. |37| However, since silicon has two atoms in its unit cell, 
successful attachment-detachment processes at the kink sites of steps must involve two 
atoms. Single atom processes from energetically favorable step edge configurations lead to 
unfavorable ones and are expected to have very large barriers. Further, the probability of 
simultaneous attachment (detachment) of two atoms is proportional to 9 2 ((1 — e ) 2 ). The 
currents are therefore given by: 

-2Eq -2(Eq+E a ) 

. vxe^Q 2 -q(1 -6) 2 . _ v u e k ^ 6 2 - c u (l - 9) 2 

Jl ? Ju 



where an is the lattice constant along the step and Ui jU are "effective" attachment attempt 
frequencies. The vanishing of these currents in equilibrium (when 9 = 9 e ) determine the 
constants q jU . Further assuming v\ = v u = v and neglecting terms of the order (9 — 9 e ) 2 
gives: 

-2E a -2(E a +E a ) 

9u9 e k b T 9u9 p H T 

a\\(l-9 e ) a\\(l-9 e ) 
from which Pq and P± can be identified. Hence 



Po + Pi 



2ua 



2 



9 



2 

e 

9~ P 



-2(Eq + E S ) 



(e k b T +e k b T ) (12) 



^ = e^ (13) 
Pi 
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Relationship between the step kinetic coefficient and the step mobility |p8| 



To relate the step mobility T and the step kinetic coefficient (3 S = Po + /3i, the expression 
for the step velocity = a 2 (j u — ji)) must be obtained in terms of (3 S . From Chernov's 
defining equations (11): 

^ = lW{X + ) + mX-) - fc9 e ) (14) 

This equation needs to be extended to account for the two dimensional character of the step: 
The currents towards the step edge are proportional to the deviation of the chemical potential 
\i from the equilibrium chemical potential near the step edge \i s e . These relationships may 
be assumed to be the defining equations for the step mobilities Ti and r : 

Assuming the existence of a coarse-grained free energy functional T of the step configuration 



{x(y)} } jj, e s is given by Mullins (43| to be: 



a 2 SFjxjy)) 

^ = ^ + W^x— (16) 

where [i e is the equilibrium chemical potential on the terrace far away from the step edge. 
Further adopting his model of a dissolving substance for the adatoms, and in an approxima- 
tion in which the deviation of the density 9 from the equilibrium density 9 e (the density far 
away from the step edge) is much smaller than 9 e , the chemical potential and the density 



can be related by p3[ : 



fx(6) -» e = hT 6 —^ (17) 

Using equations 16 and 17 it can be seen that the equilibrium density beside the step edge 
is different from 9 e . []4~Q-4~^| Further using equation 15, the step velocity can be given by: 



l4(^) + > ( ,-)-^( 1+ ^^)) 
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where T = To + ri. In the absence of an external force on the adatoms and when attachment- 
detachment of atoms from the terrace onto the step edge is rate limiting in relation to terrace 
diffusion it is expected that 6(x + ) = 8(x~) = 9 e . Then, equation 18 reduces to the Langevin 
equation (without the noise term () p9| , |3^] : 

dx T 8F(x{y)) 



dt k h T 5x 



+ C(y,t) 



Within the small slope approximation (|| << 1), ^^f^ oc . The spatial average (over 
y, denoted by <>) of this term vanishes due to periodic boundary conditions. Identifying 
< x > with X, assuming that the density is independent of y and approximating < 6(x) >= 
9(X), equation 18 (on spatial averaging) gives: 

f (19) 



Equations 14 and 19 agree if 



A),i = ^r P' = ^2 ( 2 °) 
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FIGURES 

FIG. 1. One bilayer of the Si(lll) surface consisting of the upper monolayer (grey) and lower 
monolayer (black). The figure shows the threefold and reflection symmetry of this surface: Steps 
running along directions with equal 8 are identical. (Borrowed from reference [21]) 

FIG. 2. The Si(lll) surface (a) (the upper monolayer is shown in grey and the lower monolayer 
in black), the [211] step (d) and the [TT2] step (g) (for the step configurations the upper bilayer is 
shown with larger atoms as compared to the lower bilayer). (b), (e) & (h) show the corresponding 
adatom binding energy V(x,y). (c), (f) & (i) show the corresponding binding energy derived from 
the local geometry Vi g (x,y). Minima, saddle points and maxima are marked (labeled) by +(m), 
*(s) and x(M) respectively with the figure in parenthesis being their corresponding value in eV. In 
(b) Sz (c) the contours are 0.1 eV apart. In (e), (f), (h) & (i) they are 0.2 eV apart. In (e) & (h) 
contours of V >-2.2 eV and in (f) & (i) those of V\ g >-1.9 eV are marked with dashed lines. The 
contour plots suggest a strong correlation between V and Vi g . The diffusion barrier on the surface 
is determined by ml & si in (b) & (c). The Schwoebel barrier is determined by s3 in (e) & (f) 
whereas it is determined by s5 in (h) and s6 in (i). There is a one-to-one correspondence between 
the barrier determining features in V(x, y) and Vi g (x,y) for the Si(lll) surface and the [211] step 
suggesting that the Schwoebel barrier on the [211] step is robust feature. 

FIG. 3. Plots of the adatom potential energy derived from the local geometry of fixed atoms 
Vi g vs the actual adatom potential energy V. The straight line fits show the strong correlation 
between Vi g and V. The rough configuration independence of the relationship suggests that it is a 
characteristic of the (111) surface. 
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FIG. 4. Inset (a) shows schematically the variation of r _1 = k + — pk_ with F. The solution 
to r" 1 = (of interest) is F r which exists even in the limit wX" 1 = whereas the other two do 
not. Inset (b) shows schematically the variation of F r with E s for p = 1.0 and p = 0.2. The upper 
bound on E s (E%) is determined from the condition F(E S > E^,p = 1.0) < -F(E S >E^,p = 0.2). 
It can be seen that F(E%, p = 0.2) is the lower bound on F (F e ). F = F(E S = 0,p = 0.2). For the 
curves shown, w = 0.15pm, 0.5pm < A < 0.5 x 10 5 pm (7 values equally spaced on the logarithmic 
scale) with the range of t depending on A: 7 - 7 ^ 10 ' 5 < t < 7 - 7 ^l° 3 x 10 3 - 1 (A > 7pm). The upper 
(lower) limits of t are marked by + (x). For A < 7pm the upper limit of t is forced to be smaller 
than that estimated since the solution F r ceases to exist. For A < lpm, F r does not exist for even 
the smallest t. This shows that \ m in = lpm. For wX' 1 — > 0, E^ = 0.35/c5T=0.05 eV (for the 
experiment in ref. 5) - significantly lesser than that calculated here for the [211] step (0.61±0.07 
eV). This suggests that [TT2] steps were observed. 
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Fig 2d. "Calculated Schwoebel barriers ..." by S. Kodiyalam et al. 
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Fig 2g. "Calculated Schwoebel barriers ..." by S. Kodiyalam et al. 
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Fig 2i. "Calculated Schwoebel Barriers " by S. Kodiyalam et al 
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